"""
Title: RunDoctorSim-v2.jl
Author: WB Macleod
Created: 2019-08-28
Last Edit: 2019-10-13
Description: Simulation Parmeters for Currie-Macleod paper on Doctor Decision
Making.  This is the last version for the paper. It contains code for
both 11 drugs and polypharmacy to produce a single df for ggplot2
 
This file has the parameters and does the main runs.
To use the file:

1. Make sure that one has the following files in the same directory:
     A. "DrugEffects.csv" - the mean and variance of drugs from RCT evidence.
     B. "RunDoctorSim-v2.jl" - this file
     C. "DoctorSim.jl" - Functions used for the simulation.
2. In terminal do (code has been tested on julia 1.2):
     julia RunSim-vs.jl
3. This produces two files:
    A. CM-Sim-v2.csv - simulation data
    B. CM-Sim-v2.md  - text file with parameter files.

Near then end of the file there is code that can be un-commented that will
produce time stamped files to allows experimentation and storing of data.

Copyright 2019 by W. B. MacLeod, wbmacleod@gmail.com
This file is part of "CM-DoctorSim".

    "CM-DoctorSim" is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    "CM-DoctorSim" is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with "CM-DoctorSim".  If not, see <https://www.gnu.org/licenses/>.

"""

global Gr = 0; #Run number
global Gm = 0; #Current date (Month) of treatment global Gj = 0; #Doctor Treating patient in simulation
global Gi = 0; #Patient being treated 

## Run the helper functions;
include("DoctorSim-v2.jl")
df = DataFrame(
    Run = Int64[],       # Run number - each run does all patients and docs. New seeds on each run.
    Month = Int64[],
    Did = Int64[],  #Doc index from list
    Dstyle = String[],    #long run or short run
    Dqual = String[],
    Entropy = Float64[],
    Utility = Float64[],
    Popt = Float64[], #Fraction of patients getting the best drug
    Poly = Int[],  #1 if polypharmacy is allowed.
    PolyOpt = Float64[] #Fraction of patients where polypharmacy is optimal
)
# Doctor Styles - fix reservation proability
d_style = Dict("Low-RP" => 0.5, "High-RP" => 0.95)
# This is the ability - precision of observation
d_qual = Dict("High" => 1, "Medium" => 0.10, "Low" => 0.01)
## In polypharmacy e = ppcom(e_A + e_B) + pperr * N(0,1)
ppcom = 0.6
pperr = 0.01
## There is an unobserved fixed effect for each patient
## var = pcorr * min(var), where pcorr in (0,1)
pcorr = 0.9
## Parameters for the simulation. 
NP = 300 #number of patients
RUNS = 10 #number of runs
NMONTHS = 36 #Months in a Run
"""
Next creates a list of:
    1. patients (P)
    2. Doctors to treat them
    3. Record of each visit that is updated to show knew beliefs

The idea is to have each doctor treat the same group of randomly
drawn patients.  Since we start at the same place each time.
In terms of the seed, lets reseed the random number generator
at the beginning for each doctor.  Hence, two doctors with the
same characteristics will get the same outcomes.
"""

"""
The function Treat!(D,P,PD) takes doctor, patient and record,
updates record and chooses next treatment.
Next suppose  idefine function that takes a single doctor:
Observes patient, updates beliefs, chooses drug.
Data recorded to global database
"""
#==========================================================
# Simulation with no polypharmacy
# Drug 1 is placebo
# Set beliefs and start simulations
==========================================================#
global seed=2019
mu0, Sig0 = DataIn(pcorr);
ndrugs = length(mu0)
for Sim in 1:RUNS
#Sim =1
    global Gr = Sim;
    global seed = nextprime(seed+10 ) ;  #Change seed for each run 
    Random.seed!(seed)
    P = P_list(NP, mu0, Sig0); #New set of patients

#
# Doctor Loop - Doctors face Same patients on each loop
# reset seed so their diagnostic noise is identical, so
# variation is due to Doctor effects only
    Docs = D_list(d_qual, d_style)
    for j in 1:6
        global Gj = j;
        Random.seed!(nextprime(seed + 5))
        PD = PD_list(P, mu0, Sig0); #intialize belief re Patients.
        for mth in 1:NMONTHS
            global Gm = mth;
            Monthly!(Docs[j], P, PD, df);
            print("Sim Progress with NO PolyPharmacy:\n\t ", "Run: ", Sim, " Doctor: ", j, " Month: ", mth, "\n")
        end
    end
end
#==========================================================
# Simulation with polypharmacy
# Drug 1 is placebo
# Set beliefs and start simulations
==========================================================#
global seed = 2019
mu0, Sig0, index = DataInPoly(ppcom, pperr, pcorr);
npoly = length(mu0) #number of drugs with polypharamacy
# Create an index matrix to compute entropy for single drugs
ind = zeros(ndrugs, npoly)
for k in 1:npoly
    for i in index[k]
        ind[i,k] = 1
    end
end

for Sim in 1:RUNS
#Sim =1
global Gr = Sim;
    global seed = nextprime(seed+10 ) ;  #Change seed for each run 
    Random.seed!(seed)
    P = P_list(NP, mu0, Sig0, ndrugs); #New set of patients

#
# Doctor Loop - Doctors face Same patients on each loop
# reset seed so their diagnostic noise is identical, so
# variation is due to Doctor effects only
#
    Docs = D_list(d_qual, d_style)
    for j in 1:6
        global Gj = j;
        Random.seed!(nextprime(seed + 5))
        PD = PD_list(P, mu0, Sig0); #intialize belief re Patients.
        for mth in 1:NMONTHS
            global Gm = mth;
            #Monthy computes stats and saves 1 at end indicates polypharmacy, ind is to adjust entropy
            Monthly!(Docs[j], P, PD, df, 1, ind);
            println("Sim Progress with PolyPharmacy:\n\t ", "Run: ", Sim, " Doctor: ", j, " Month: ", mth, "\n")
        end
    end
end

#Save the results
Time = Dates.format(now(),"Y-m-d-HH-SS");
#File = "CM-Sim-v2-" * Time  #uncomment this line and comment following line to have time stamped output.
File = "CM-Sim-v2"
CSV.write(File * ".csv", df);
## Write out parameter values for Run.
## Use markdown for headers etc

open(File * ".md", "a") do f
    println(f, "# Parameters Values for Simulation at Time: ", Time)
    println(f, "## Doctor Parameters")
    println(f, "### Doctor's Reservation Probability:")
    for (k,v) in d_style
        println(f, k, ": ", v)
    end
    println(" ")
    println(f, "### Doctor's Skill:")
    for (k, v) in d_qual
        println(f, k, ": ", v)
    end
    println(f, "## Patient Parameters")
    println(f, "Degree of Complementarity in PolyPharmacy: ", ppcom)
    println(f, "Idiosyncratic Variance in PolyPharmacy: ", pperr)
    println(f, "Number of Patients per Doctor: ", NP)
    println(f, "## Simulation Parameters")
    println(f, "Months of Treatment: ", NMONTHS)
    println(f, "Number of Simulation Runs: ", RUNS)
end
